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Abstract: 

In this article, we present a galactic gravitational model of three degrees of freedom (3D), in order to study and 
reveal the character of the orbits of the stars, in a binary stellar system composed of a primary quiet or active 
galaxy and a small satellite companion galaxy. Our main dynamical analysis will be focused on the behavior of 
the primary galaxy. We investigate in detail the regular or chaotic nature of motion, in two different cases: (i) 
the time-independent model in both 2D and 3D dynamical systems and (ii) the time-evolving 3D model. For the 
description of the structure of the 2D system, we use the classical method of the Poincare (x, p x ), y = 0, p y > 
phase plane. In order to study the structure of the phase space of the 3D system, we take sections in the plane y = 
of the 3D orbits, whose initial conditions differ from the plane parent periodic orbits, only by the z component. The 
set of the four-dimensional points in the (x, p x , z, p z ) phase space is projected on the (z, p z ) plane. The maximum 
Lyapunov characteristic exponent is used in order to make an estimation of the chaoticity of our galactic system, 
in both 2D and 3D dynamical models. Our numerical calculations indicate that the percentage of the chaotic orbits 
increases when the primary galaxy has a dense and massive nucleus. The presence of the dense galactic core also 
increases the stellar velocities near the center of the galaxy. Moreover, for small values of the distance R between 
the two bodies, low-energy stars display chaotic motion, near the central region of the galaxy, while for larger 
values of the distance R, the motion in active galaxies is entirely regular for low-energy stars. Our simulations 
suggest that in galaxies with a satellite companion, the chaotic nature of motion is not only a result of the galactic 
interaction between the primary galaxy and its companion, but also a result caused by the presence of the dense 
nucleus in the core of the primary galaxy. Theoretical arguments are presented in order to support and interpret 
the numerically derived outcomes. Furthermore, we follow the 3D evolution of the primary galaxy, when mass 
is transported adiabatically from the disk to the nucleus. Our numerical results are in satisfactory agreement with 
observational data obtained from the M51-type binary stellar systems. A comparison between the present research 
and similar and earlier work is also made. 
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1 Introduction 

It is well known that galaxies are not observed as isolated 
entities, but often appear in binary, triple, or multiple stellar 
systems. On the contrary, just as there are clusters of stars, 
there are also clusters of galaxies. This leads to gravitational 
galactic interactions. The Milky Way-Magellanic Clouds 
system is a very indicative example of galactic interaction 
(see Gardiner, Sawa, & Fujimoto 1994; Lin, Jones, & Kle- 
mola 1995; Weinberg 1995, 1998; Van der Marel 2001; Put- 
man et al. 2003; Connors et al. 2004; Bekki & Chiba 2005). 
Another interesting stellar system of interacting galaxies is 
the Andromeda galaxy M3 1 , with its two smaller companion 
galaxies M32 and NGC 205. The gravitational effects of the 
M32 galaxy on the spiral structure of M31 were studied by 
Byrd (1978). This interesting triple stellar system was also 
investigated in an earlier work, using a self-consistent com- 
puter simulation code with revealing results (see Vozikis & 
Caranicolas 1994). 

One of the main reasons for the transition from ordered 
to chaotic motion is the presence of massive objects in the 



central regions of the galaxies. Stars reaching the galactic 
center on highly eccentric redial orbits can be scattered out 
of the galactic plane displaying chaotic motion (see Sellwood 
& Moore 1999). Furthermore, a central mass concentration 
can strongly perturb the stellar box orbits in elliptical galax- 
ies which become chaotic (Merritt 1996). Observational in- 
dications suggest the presence of a strong central mass con- 
centration by a very sharply rising rotation curve. A second 
reason for chaotic behavior in the galactic motion is the pres- 
ence of strong external perturbations (see Zotos 2011c). In- 
creasing the magnitude of perturbations weakens or destroys 
the stability of orbits, thereby increasing the percentage of 
the stochasticity in the stellar system. Resonances are also 
responsible for chaotic motion (see Contopoulos & Grosb0l 
1986; Cincotta, Giordano, & Perez 2006). Of course, res- 
onances are responsible not only for the presence of chaos 
in galactic stellar systems but also for the chaotic motion in 
the solar system (see Wisdom 1987; Henrard & Caranicolas 
1990). The reader can find interesting information about the 
chaotic motion in galaxies and its connection with observa- 
tional data in Grosb0l (2002). 
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An extended and diffuse stellar halo envelops the Milky 
Way. Although only an extremely small fraction of the stars 
in the Solar neighborhood belong to this halo, they can be 
easily recognized by their extreme kinematics and metallici- 
ties. Stellar populations with these properties can now be fol- 
lowed to distances in excess of 100 kpc using luminous trac- 
ers such as RR Lyraes, blue horizontal branch stars, metal- 
poor giants, and globular clusters (e.g. Oort 1926; Baade 
1944; Eggen, Lynden-Bell, & Sandage 1962; Searle & Zinn 
1978; Vivas & Zinn 2006; Morrison et al. 2009). 

Low surface-brightness features seen in projection around 
other galaxies aid in the interpretation of the Milky Way's 
stellar halo and vice versa. Diffuse concentric 'shells' of stars 
on 100 kpc scales around otherwise regular elliptical galax- 
ies have been attributed to accretion events (e.g. Schweizer 
1980; Quinn 1984). Recent surveys of M31 (e.g. Ferguson 
et al. 2002; Kalirai et al. 2006; McConnachie et al. 2009) 
have revealed an extensive halo (to ~ 150 kpc) also display- 
ing abundant substructure. The surroundings of other nearby 
Milky Way analogues are now being targeted by observa- 
tions using resolved star counts to reach very low effective 
surface brightness limits, although as yet no systematic sur- 
vey has been carried out to sufficient depth (e.g. Zibetti & 
Ferguson 2004; McConnachie et al. 2006; de Jong, Radburn- 
Smith, & Sick 2008; Barker et al. 2009; Ibata, Mouhcine, & 
Rejkuba 2009). A handful of deep observations beyond the 
Local Group suggest that stellar halos are ubiquitous and di- 
verse (e.g. Sackett et al. 1994; Shang et al. 1998; Malin & 
Hadley 1999; Martmez-Delgado et al. 2009; Faundez-Abans 
et al. 2009). 

Stellar halos formed from the debris of disrupted satel- 
lites are a natural byproduct of hierarchical galaxy forma- 
tion in the Lamda Cold Dark Matter (ACDM) cosmology. 
The entire assembly history of a galaxy may be encoded 
in the kinematics, metallicities, ages, and spatial distribu- 
tions of its halo stars. Even though these stars constitute a 
very small fraction of the total stellar mass, the prospects are 
good for recovering such information from the halos of the 
Milky Way, M31, and even galaxies beyond the Local Group 
(e.g. Johnston, Hernquist, & Bolte 1996; Helmi & White 
1999). In this context, theoretical models can provide useful 
'blueprints' for interpreting the great diversity of stellar ha- 
los and their various sub-components, and for relating these 
components to fundamental properties of galaxy formation 
models. Alongside idealized models of tidal disruption, ab 
initio stellar halo simulations in realistic cosmological set- 
tings are essential for direct comparison with observational 
data. 

Therefore, it seems very challenging to construct a galac- 
tic gravitational model of three degrees of freedom (3D), in 
order to study and reveal the character of the orbits of the 
stars, in a binary stellar system composed of a primary quiet 
or active galaxy and a small satellite companion galaxy. The 
aim of the present paper is to study the character (regular 
or chaotic) of motion in the primary galaxy (hereafter the 
galaxy) under the perturbation caused by a satellite galaxy 
and to connect the degree of chaos with parameters such as 
the mass of the nucleus of the galaxy or the distance between 
the two galaxies. We shall also study the behavior of the 
velocities near the central region and try to connect this be- 
havior with the scale length of the nucleus. Furthermore, we 
shall deal with the case where adiabatic mass transfer takes 
place in the galaxy, that is, the galaxy evolves with time. In 



this initial research, only the case where the two bodies move 
in circular orbits about the center of mass of the system is 
considered, for simplicity. We plan to study the case of the 
satellite in inclined, elliptic orbits in a future work. 

In a recent paper (Zotos 2012b), we investigated in de- 
tail and revealed the regular or chaotic character of motion 
in a galactic gravitational model of three degrees of freedom 
describing a binary quasar system. According to the model, 
two quasars are hosted in a pair of interacting disk galaxies. 
Here we must emphasize that in the present research the sec- 
ond body of the binary stellar system, that is, the satellite 
companion galaxy, is considered and treated as a point mass. 
On the contrary, in Zotos (2012b), we deal with a much more 
complicated binary system, since none of the two interacting 
disk galaxies can be treated as point mass due to the fact that 
both of them are hosting quasars in their cores. Thus, we may 
say that Zotos (2012b) is a reasonable extension and gener- 
alization of the present paper. Taking into account that in 
both papers we study the properties of motion of the stars in 
similar, but in no case identical binary systems, it is natural 
and well expected to use the same philosophy regarding the 
basic setup of the gravitational models, the handling of the 
numerical calculations, and also the same standard methods 
of galactic dynamics and celestial mechanics. 

The present article is organized as follows. In Section 
2, we present our gravitational dynamical model, which de- 
scribes the motion in a binary stellar system. In Section 3, 
we provide an analysis of the two degrees of freedom (2D) 
system, considering orbits in the galactic plane (z = 0). In 
Section 4, we study the character of motion in the 3D sys- 
tem, using different kinds of dynamical methods. In Section 
5, we present an interesting theoretical analysis in order to 
support the numerically obtained results. In Section 6, we 
use a 3D time-dependent model in order to follow the evo- 
lution of orbits as mass is transported adiabatically from the 
disk to the nucleus of the primary galaxy. We close this pa- 
per with Section 7, where the conclusion is presented and a 
comparison between the present theoretical results with ob- 
servational data is made. In the same section, we include 
statements concerning the astrophysical relevance of this re- 
search work. 

2 Presentation of the Dynamical Model 

Our gravitational model describes the dynamical properties 
of a binary stellar system which is composed of a primary 
galaxy and a small satellite companion galaxy. The primary 
galaxy (hereafter the galaxy) consists of a disk and a spheri- 
cal nucleus and is described by the potential 

V P (r,z) = V i (r,z)+V n (r,z) 

GMd 

^b 2 + r 2 + (a + tJWTz 1 ) 2 



where r 2 = x 2 + y 2 , M A and M n is the mass of the disk and 
the nucleus of the primary galaxy respectively, a is the disk's 
scale length, h is the disk's scale height, b is the core radius 
of the disk halo, while c„ is the scale length of the nucleus. 
The disk of the primary galaxy is represented by the well- 
known Miyamoto-Nagai model (Miyamoto & Nagai 1975), 
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with an additional core radius parameter (see Carlberg & In- 
nanen 1987). The additional core radius is used so that the 
rotation curve displays high values at large distances from 
the galactic center, justifying in a way the presence of dark 
matter. The Plummer sphere we choose to describe the nu- 
cleus has been used many times in the past, in order to study 
the effects of the introduction of a central mass component 
in the core of a galaxy (see Hasan & Norman 1990; Hasan et 
al. 1993). The satellite galaxy is described by a point-mass 
potential 

V«. d =™ (2) 
d 

where d 2 = x 2 + y 2 + z 2 , while M s is the mass of the satellite 
galaxy. Therefore, the total gravitational potential describ- 
ing the motion in our binary system is Vj(r, z) = Vp(r, z) + 

v s (nz). 

In our study, we shall use the theory of the circular re- 
stricted three-body problem. The two bodies move in circular 
orbits in an inertial frame OXYZ, with the origin at the cen- 
ter of mass of the system, with a constant angular frequency 
n p > 0, given by Kepler's third law 



GM t 



(3) 



where M, = M A + M n + M s is the total mass of the binary 
system, while R is the distance between the centers of the 
two galaxies. A clockwise, rotating frame Oxyz, is used with 
the axis Oz coinciding with the axis OZ and the axis Ox 
coinciding with the straight line joining the two bodies. In 
this frame, which rotates with angular frequency Cl p , the two 
galactic centers have fixed positions C\(x,y, z) = (xi,0, 0) 
and Ci{x,y, z) = (x 2 , 0, 0) , respectively. The total potential 
which is responsible for the motion of a star in the dynamical 
system of this binary system is 

<D t (x,y,z) = 3>p(x,y,z) + ®s(x,y,z) + <J> rot (x,y), (4) 



where 

<S>t>(x,y,z) 



GM,\ 



yjb 2 + r 2 al + (a + WVz 2 ) 2 
GM„ 



®s(x,y,z) 



GM S 



<S>rot(x,y) 



M,' 



,M S 
R 2 —R, 



(5) 



and 



(x 



xxf+y 2 , 
x 2 ) 2 +y 2 , 



+ z 



(6) 
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Mt' 



(7) 



The angular frequency Q. p is calculated as follows. The 
two bodies circulate around their common mass center of the 
system with angular frequencies Q. pl and Sl p2 given by 



i_(-dV s (r) 

xi 



dr 



n 



l'2 



]_ /dV P (r) 

Xt 



dr 



(8) 



As the primary galaxy cannot be considered as a mass point, 
the two angular frequencies are not equal, in general. How- 
ever, this issue can easily be resolved. The angular frequen- 
cies of the two bodies can become equal under reasonable 
assumptions. This is justified if the final set of the param- 
eters has physical meaning and represents satisfactorily the 
dynamical system. The equation Q pl = Q p2 leads to a four- 
fold infinity of solutions in the four unknowns (a, b, h, c n ). If 
one chooses proper values (representing the dynamical sys- 
tem), say, for the three parameters (b, h, c„), then the equation 
Q„i = fl P 2 gives only two values for the parameter a. One 
value is positive, while the other is negative and is rejected. 
The author would like to make it clear that after choosing 
properly the parameters, the deviation between the two angu- 
lar frequencies is negligible, so that v = [\Q p i - SlpzU /Q p i 
is of the order of 1(T 8 or even smaller and £ = \Q p i - Cl p \ 
or £ = |£2 p2 - Q p \ is of the order of 10~ 6 . Therefore, we 
consider the two angular frequencies almost equal, that is 
Q. P i = Q. P 2 = ii p . Moreover, the treatment of large bodies 
with spherical symmetry as mass points is very common in 
celestial modeling. This method is quite familiar to those 
measuring the masses of disk galaxies. 

The above setup in the rotating frame has been applied 
successfully in several previous papers studying the nature 
of binary stellar systems (see Caranicolas & Papadopoulos 
2009; Caranicolas & Zotos 2009a; Caranicolas & Innanen 
2009; Zotos 2012b). Thus, we may conclude that the dy- 
namical model described in Equations (4) and (5) is a good 
and also realistic model for a pair of interacting galaxies at 
the given mass ratio. Furthermore, in a recent paper (Zo- 
tos 2012b), we used the same treatment and setup, in order 
to study the dynamical properties in a binary system of two 
interacting galaxies. This could be considered a reasonable 
extension of the present study, where the second galaxy here 
is assumed to be only a satellite companion. The fact that 
the primary galaxy is sufficiently apart from its satellite al- 
lows us to assume that the tidal phenomena are very small 
and therefore negligible. In particular, our numerical calcula- 
tions indicate that the tidal phenomena are significant enough 
when the distance between the centers of the two galaxies is 
R < 0.95. In the present study, we assume that the range of 
this distance is 1.5 < R < 3 and, therefore, the tidal forces 
can be neglected. 

In this rotating frame, the equations of motion are 



2Q. p y, 



dx 

-? + 2<V, 

dy 

d<P, 



(9) 



where the dot indicates derivative with respect to the time. 
The only integral of motion for the system of differential 
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Figure 2: Contours of the projections of the isopotential 
curves O t (x, y, z) = £/ on the (x, y) plane. The five Lagrange 
equilibrium points are indicated as Li, Li, L3, L4, and L5, 
while Ci and C2 are the centers of the two galaxies at a dis- 
tance of R = 1.5. 



equations (9) is the well-known Jacobi integral given by the 
equation 

J=\ [pI + p\ + Pi) + 3>t(*,.y, z) = Ej, (10) 

where p x , p y , and p z are the momenta per unit mass conjugate 
to x, y, and z, while Ej is the numerical value of the Jacobi 
integral. 

All numerical outcomes of the present work are based 
on the numerical integration of the equations of motion (9), 
which was made using a Bulirsh-Stoer routine in fortran77, 
with double precision in all subroutines (see Numerical Recipes 
in FORTRAN, 2nd edn, in Press et al. 1992). The accu- 
racy of the calculations was checked by the consistency of 
the Jacobi integral (10), which was conserved up to the 12th 
significant figure. 

In this article, we shall use a system of galactic units 
where the unit of length is 20 kpc, the unit of mass is 1.8 x 
10" Mq, and the unit of time is 0.99 x 10 s yr. The velocity 
unit is 197 kms~', while G is equal to unity (see Vozikis & 
Caranicolas 1992). In these galactic units, we use the val- 
ues a = 0.15, b = 0.368, h = 0.00625, c n = 0.0125, and 
M s = 0.2. The values of the above quantities of the dynam- 
ical system remain constant during this research, while the 
values of M d , M„, and R are treated as parameters. The above 
numerical values of the constant dynamical quantities of the 
system secure positive density everywhere and free of singu- 
larities. 

Figures l(a-c) show the contours of the projections of 
the isopotential curves O t (x, y, z) = Ej on the (x, y), (x, z) and 
(y, z) primary planes, respectively. The values of the param- 
eters are M A = 2, M„ = 0.08, R = 1.5, and Sl p = 0.821922. 
Lighter colors indicate higher values of Ej. Figure 2 shows 
the contours of the isopotential curves 0,(x, y, z) = Ej, on 
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Figure 3: The rotation curve of the primary galaxy is shown 
as the black line. The red line is the contribution from the 
spherical nucleus, while the blue line is the contribution from 
the disk-halo potential. (Colors available only in the elec- 
tronic version of the article). 



the (x,y) plane, with some additional details regarding the 
stability and the structure of the dynamical system. More- 
over, Li, L2, Li, L4, and L5 are the five Lagrange equilibrium 
points, while C\ and Ci are the centers of the two galaxies at 
a distance of R = 1.5. At these equilibrium points, we have 

dO t (x,y,z) = dO,(x,y,z) = d<t> t (x,y,z) = Q 
dx dy dz 

Note that Lj, L 2 , and L 3 are the unstable saddle equilibrium 
points, while L 4 and L 5 are the triangular points (see Binney 
& Tremaine 2008). 

The rotation curve of the primary galaxy when M d = 
2 and M n = 0.08 is shown as the black line in Figure 3. 
In the same plot, the red line is the contribution from the 
spherical nucleus, while the blue line is the contribution from 
the disk-halo component. We observe that at small distances 
from the primary galactic center r < 0. 1 length units, which 
equals 2 kpc, dominates the contribution from the spherical 
nucleus, while at larger distances r > 2 kpc, the disk-halo 
contribution is the dominant factor. 

It would be very illuminating to compute the mass den- 
sity distribution p(x, y, z) of our binary system, derived from 
the total potential V T (r, z), using the Poisson's equation 

P(x,y,z) = -^-V 2 V T (x,y,z) 
An 

1 (d 2 V T 8 2 V T d 2 V T \ 
An \ dx 1 + dy 2 + dz 2 )' 

Figures 4(a-c) show the projections of the isodensity curves 
p(x, y, z) = const on the (x, y), (x, z), and (y, z) primary planes, 
respectively, when M A = 2, M n = 0.08, R = 1.5, and Q. p = 
0.821922. In the first two plots, we can observe the small 
variations on density distribution, caused by the mass of the 
satellite galaxy. Figures 5(a-c) depict the 3D plots of the val- 
ues of the total mass density p(x, y, z) on the same primary 
planes. Once again, we can distinguish the spherical nucleus 
at the central region of the primary galaxy, the well-formed 
disk structure, and also the abrupt peak of the mass density 
at the position of the satellite galaxy (x 1.36). These 3D 
plots ensure that the mass density of our binary galactic sys- 
tem has positive values everywhere. Going one step further, 



Figure 1: (a-c) Contours of the projections of the isopotential curves 0,(x, v, z) = £/ on the (x,y), (x, z), and (y,z) planes. 
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Figure 4: (a-c) Contours of the projections of the isodensity curves p(x, y, z) = const on the (x, y), (x, z), and (y, z) planes. 




Figure 5: (a-c) 3D plots of the distribution of the total mass density p(x,y,z) on the (x,y), (x,z), and (y,z) planes. 




Figure 6: (a-c) The evolution of the total mass density along the x, y, and z axes. 



we present in Figures 6(a-c) the evolution of the mass den- that in all cases the mass density obtains high values near 
sity along the x, y, and z axes, respectively. One can observe the two galactic centers, while it reduces and tends asymp- 
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totically to zero, as the distance from the centers of the two 
galaxies increases. The values of the mass density p on the 
vertical axis are in M Q pc~ 3 . 

3 Structure of the 2D Hamiltonian 

System 

In this section, we shall investigate the properties of motion 
in the Hamiltonian system of two degrees of freedom. This 
can be derived from Equation (10) if we set z = p z = 0. 
Then, the corresponding Hamiltonian is 

h = \ (pI + p;) + ® t (x,y) = E J2 , (13) 

where E J2 is the numerical value of J 2 . As the dynamical sys- 
tem is now two-dimensional, we can use the classical, qual- 
itative method of plotting the successive intersections of the 
2D orbits, using the (x, p x ), y = 0, p y > Poincare surface 
of section (PSS), in order to determine the regular or chaotic 
character of motion. This method has been extensively ap- 
plied to Hamiltonian systems with two degrees of freedom, 
as in these systems the PSS is a two-dimensional plane. The 
results obtained from the study of the 2D system will be used 
to help us understand the structure of the more complicated 
phase space of the 3D system, which will be presented in the 
following section. 

Figures 7(a-b) show the (x, p x ), y = 0, p y > Poincare 
phase plane, for the motion of a star in the primary galaxy, 
which obtained by numerical integration of the equations of 
motion (9). In this case, the distance between the centers of 
the two galaxies is R = 1.5 and Q p = 0.821922. The value 
of the Jacobi integral is Ej2 = -2.6. This particular value of 
the Jacobi integral will remain constant for all our numerical 
experiments. In Figure 7(a) we have M A = 2 and M„ = 0.08, 
that is, the primary galaxy has a dense and massive nucleus. 
We observe a large chaotic area, while the regular region is 
mainly confined around the stable retrograde periodic point. 
There are also some small islands of invariant curves em- 
bedded in the chaotic sea, corresponding to secondary reso- 
nant orbits. Note that near the center of the primary galaxy, 
the velocity obtains high values of the order of 900 kms~'. 
This results from the presence of the dense nucleus in the 
core of the galaxy and it is characteristic of strong galactic 
activity. The outermost black solid line defines the zero ve- 
locity curve. Figure 7(b) is similar to Figure 7(a), but when 
M d = 2.08 and M u = 0. As the total mass of the system 
is conserved, this means that now the disk occupies all the 
available mass of the primary galaxy. In this case, the pattern 
has two main differences from the pattern shown in Figure 
7(a). The first difference is that the chaotic region is smaller 
and the second is that the velocity near the galactic core is 
smaller, about 470 kms~', than that of Figure 7(a). More- 
over, some secondary resonances appear in the area of the 
retrograde orbits. Some small sticky regions are also present 
in both cases (active or quiet galaxy). Therefore, one can say 
that the structure of chaos in galaxies is a result not only of 
galactic interaction but also of the nuclear galactic activity. In 
other words, the presence of the dense and massive nucleus 
in the galactic core of the primary galaxy affects drastically 
the character of motion of the stars. 

Figure 8(a) is similar to Figure 7(a) but when the dis- 
tance between the centers of two the galaxies is R = 2 and 
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Figure 1 1 : A plot of the area A% of the (x, p x ) phase plane 
covered by chaotic orbits as a function of the distance R, 
when the primary galaxy is active or quiet. 



Q.,, = 0.533854. The majority of the orbits around the direct 
and retrograde periodic point are ordered. There is a signifi- 
cant chaotic area mainly near the central region of the galaxy. 
A high-velocity value is once again observed near the galac- 
tic core of the primary galaxy, while some secondary reso- 
nances are also present. Figure 8(b) is similar to Figure 8(a), 
but when the galaxy is quiet, that is, when M n = 0. Here, 
the secondary resonances look more prominent, while only 
a small and confined chaotic layer around the direct periodic 
point is present. Furthermore, as the primary galaxy is quiet, 
the velocity near the center of the galactic core is smaller 
than that observed in Figure 8(a). Figures 9(a-b) are similar 
to Figures 8(a-b), but when the distance between the centers 
of the two galaxies is R = 2.5 and fl p = 0.381995. In Figure 
9(a), we observe that there is a relatively small chaotic layer 
confined mainly near the central region of the galaxy, while 
the rest of the phase plane is covered by invariant curves cor- 
responding to regular orbits circulating around the stable di- 
rect and retrograde periodic point. Figure 9(b) is similar to 
Figure 9(a) when the galaxy is quiet. All orbits seem to be 
regular. Chaotic motion was not observed and, if present, is 
negligible. Figures 10(a-b) are similar to Figures 9(a-b) but 
when the distance between the centers of the two galaxies is 
R = 3 and Q p = 0.290593. In Figure 10(a), we observe that 
the entire phase plane is covered by invariant curves corre- 
sponding to regular orbits circulating around the stable direct 
and retrograde periodic point, while there is no indication of 
chaotic motion. Figure 10(b) is similar to Figure 10(a), but 
when the galaxy is quiet. Once again, all orbits are regular 
and chaotic motion was not observed. Therefore, we may 
conclude that the primary galaxy, for large separations from 
its companion, does not show chaotic motion irrespective of 
the presence of the dense and massive nucleus in the galactic 
core. The outermost solid curve in the (x, p x ) phase planes 
shown in Figures 7-10 is the limiting curve in each case. 

Figure 1 1 shows the percentage A% on the (x, p x ) phase 
plane covered by chaotic orbits as a function of the distance 
R between the centers of the two galaxies, in two different 
cases. In the first case, the primary galaxy is active and the 
values of the parameters are M A = 2 and M n = 0.08, while 
in the second case the primary galaxy is quiet and, therefore, 
the values of the parameters are M d = 2.08 and M n = 0. The 
value of the angular frequency fl ; , is calculated for each par- 
ticular value of the distance R from Equation (3). The range 
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of the values regarding the distance between the centers of 
the two galaxies is 1.5 < R < 3. We observe in Figure 11 
that the chaotic percentage A% decreases exponentially, as 



the distance R increases in both cases (active and quiet pri- 
mary galaxy). Moreover, both fitting curves converge to zero, 
when R = 3. The lower fitting curve, which corresponds to 
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Figure 10: (a-b) Similar to Figures 7(a-b) but when R = 3 and Q p = 0.290593. 




Figure 12: A plot of the average value of the (LCE) as a 
function of the distance R, when the primary galaxy is active 
or quiet, for the 2D dynamical system. 



the case when the primary galaxy is quiet, presents a more 
rapid reduction and reaches more quickly (when R > 2.35) 
to the zero limit. A more detailed view of Figure 1 1 reveals 
that for each particular value of the distance R, the chaotic 
percentage A% is always smaller when the primary galaxy is 
quiet. We must point out that the chaotic percentage A% is 
calculated as follows: we choose 10 3 orbits with random ini- 
tial conditions (x a ,p x0 ) in eac h phase plane and then divide 
the number of those who correspond to chaotic orbits to the 
total number of the tested orbits. 

To have a better estimation about the degree of chaos in 
our dynamical system in every case, we have computed the 
average value of the Lyapunov characteristic exponent (LCE; 
see Lichtenberg & Lieberman 1992). Figure 12 shows a plot 
of the (LCE) as a function of the distance R between the cen- 
ters of the two galaxies, in two different cases. As previously, 
in the first case, the primary galaxy is active and the values of 
the parameters are Ma = 2 and M n = 0.08, while in the sec- 
ond case the primary galaxy is quiet and, therefore, the val- 
ues of the parameters are M A = 2.08 and M n = 0. The value 
of the angular frequency Q. p is calculated for each particu- 
lar value of the distance R from Equation (3). The range of 
values regarding the distance between the centers of the two 
galaxies is again 1.5 < R < 3. One can see in Figure 12 that 



(LCE) decreases almost linearly as the distance R increases 
when the primary galaxy is active. On the contrary, when 
the nucleus is absent, we observe that the corresponding fit- 
ting curve presents a more rapid reduction. Once again, as 
pointed out in Figure 1 1, for each particular value of the dis- 
tance R, the average value (LCE) is always smaller when the 
primary galaxy is quiet. Moreover, both fitting curves con- 
verge to zero when R = 3. Here we must note that it is well 
known that the value of the LCE is different in each chaotic 
component (see Saito & Ichimura 1979). As we have in all 
cases regular regions and only one unified chaotic area in 
each (x, p x ) phase plane, we calculated the average value of 
the LCE by taking 10 3 orbits with random initial conditions 
(xcp.to) in the chaotic area in each case and we integrated 
each orbit for a time interval of 10 5 time units in order to ob- 
tain a reliable value of the LCE regarding the nature of the 
orbit. Here we have to point out that all calculated LCEs cor- 
responding to chaotic orbits were different only on the fifth 
decimal point in the same chaotic region. 

In Figures 13(a-h) we present eight representative orbits 
of the 2D dynamical system. Figure 13(a) shows a quasi- 
periodic orbit circulating around the primary galaxy, with 
initial conditions x = 0.167, y = 0, and p v0 = 0, while 
the value of p v o is obtained from the Jacobi integral (13) 
for all orbits. The values of all the other parameters are as 
in Figure 7(a). Figure 13(b) shows a periodic orbit moving 
around the center of the primary galaxy, with initial condi- 
tions x'o = -0.1895, jo = 0, and p x0 = 0, while the val- 
ues of all other parameters are as in Figure 7(b). This or- 
bit is characteristic of the 1:2 resonance. In Figure 13(c), a 
figure-eight type periodic orbit, characteristic of the 1:3 res- 
onance, is shown. The initial conditions are xo = 0.5788, 
yo = 0, and p v0 = 0, while the value of all other param- 
eters are as in Figure 8(a). Figure 13(d) shows a periodic 
orbit moving around the primary galaxy with initial condi- 
tions x = 0.4295, yo = 0, and p M = 0, while the values of 
all the other parameters are as in Figure 7(b). This orbit is a 
characteristic example of the 3:3 resonance. In Figure 13(e), 
a periodic orbit with initial conditions x = 0.7445, yo = 0, 
and p x o = is presented. This orbit belongs to the family 
of the 2:3 resonant orbits. The values of all other parameters 
are as in Figure 7(a). In Figure 13(f), we see a periodic orbit 
with initial conditions x = -0.9854, y = 0, and p v0 = 0. 




Figure 13: (a-h) Eight representative regular orbits of the 2D dynamical system. The values of the initial conditions and all the 
other parameters are given in the text. 
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The values of all other parameters for this orbit are as in Fig- 
ure 7(b). Figure 13(g) shows a complicated resonant peri- 
odic orbit of higher multiplicity, circulating around the pri- 
mary galaxy, with initial conditions xo = 0.673, vo = 0, and 
p A o = 0, while the values of all the other parameters are as in 
Figure 7(a). This orbit produces a set of nine tiny islands of 
invariant curves which are embedded in the retrograde area 
of the unified chaotic domain shown in Figure 7(a). A chaotic 
orbit with initial conditions xo = -0.31, yo = 0, and p x o = 
is given in Figure 13(h). The values of all other parameters 
for this orbit are as in Figure 7(a). We observe that stars mov- 
ing in chaotic orbits pass arbitrarily around the galactic core 
without getting close enough. All orbits shown in Figures 
13(a-h) were calculated for a time period of 150 time units. 

4 Structure of the 3D Hamiltonian 

System 

In this section, we shall investigate the regular or chaotic na- 
ture of motion in the 3D Hamiltonian system described by 
Equation (10). In order to keep things simple, we shall use 
our experience gained from the study of the 2D dynamical 
system in order to obtain a clear picture regarding the prop- 
erties of motion in the 3D dynamical model. We are particu- 
larly interested in locating the initial conditions in the 3D dy- 
namical system, producing regular or chaotic orbits. A con- 
venient way to obtain this is to start from the (x, p x ) phase 
planes of the 2D system with the same value of the Jacobi 
integral as used in the 2D system and described in the previ- 
ous section. Specifically, the regular or chaotic nature of the 
3D orbits is found as follows: we choose initial conditions 
(xo, PxO, Zo), yo = Pzo - 0, such that (xq, p^) is a point on the 
phase planes of the 2D system. The points (x , p M ) lie inside 
the limiting curve 

l -p\ + * t (x) = E J2 , (14) 

which is the limiting curve containing all the invariant curves 
of the 2D system. Thus, we take Ej = E J2 . For this purpose, 
a large number of orbits (about 10 3 ) were computed with ini- 
tial conditions (x , p M , zo), where (x ,p xQ ) is a point in the 
chaotic regions of the (x, p x ) phase planes of Figures 7(a-b), 
8(a-b), 9(a-b), and 10(a-b), with all permissible values of zo 
and p z o = 0. Remember that as we are on the phase plane, we 
have vo = 0, while in all cases the value of p v o was obtained 
from the Jacobi integral (10). All tested orbits were found to 
be chaotic. Moreover, it is well known that, usually, chaotic 
domains in 2D systems lead to chaotic motion in 3D space 
when an additional degree of freedom is added, like a 'per- 
turbation'. Therefore, one may conclude that orbits with the 
initial conditions (xo,p x o) be a point in the chaotic regions 
of the 2D phase planes and, for all permissible values of zo, 
remain chaotic in the 3D system as well. 

In Figure 14, we present a plot of the average value of 
the LCE of the 3D system, as a function of the distance 
R between the centers of the two galaxies, in two different 
cases. As in the previous section, in the first case, the pri- 
mary galaxy is active and the values of the parameters are 
Mj = 2 and M n = 0.08, while in the second case the primary 
galaxy is quiet and, therefore, the values of the parameters 
are = 2.08 and M n = 0. The value of the angular fre- 
quency Q. p is calculated for each particular value of the dis- 
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Figure 14: A plot of the average value of the maximum 
(LCE) as a function of the distance R, when the primary 
galaxy is active or quiet, for the 3D dynamical system. 



tance R from Equation (3). The range of values regarding 
the distance between the centers of the two galaxies is again 
1.5 < R < 3. One can see in Figure 14 that the (LCE) de- 
creases almost linearly as the distance R increases when the 
primary galaxy is active when R < 2.25. On the contrary, 
when the nucleus is absent, we observe that the correspond- 
ing fitting curve presents a more rapid reduction. Once again, 
as pointed out in Figure 12 regarding the 2D system, for each 
particular value of the distance R, the average value (LCE) is 
always smaller when the primary galaxy is quiet. Moreover, 
both fitting curves converge earlier to zero value, than in Fig- 
ure 12. The method we used in order to compute the (LCE) 
in each case is the same as described in Figure 12. For all 
tested chaotic orbits in the 3D system, the initial value of zo 
is common and equal to 0. 1 . Here we must point out that if 
we compare the plot of the 3D system shown in Figure 14 
with that of the 2D system shown in Figure 12, we observe 
that in each case (active and quiet galaxy, respectively) the 
average values of the LCEs of the 3D system are smaller than 
those of the 2D system. Therefore, we may conclude that the 
degree of chaos in the primary galaxy is smaller when the 
dynamical system has three degrees of freedom. 

Our next step is to study the character of orbits with ini- 
tial conditions (x , p x o, Zo), yo = p z o = 0, such that (x , p x o) is 
a point in the regular regions of Figures 7(a-b), 8(a-b), 9(a- 
b), and 10(a-b). The phase space of a conservative system of 
three degrees of freedom has six dimensions, i.e. in Cartesian 
coordinates (x, y, z, x, y, z). For a given value of the Jacobi 
integral, a trajectory lies on a five-dimensional manifold. In 
this manifold, the surface of the section is four-dimensional. 
This does not allow us to visualize and interpret directly the 
structure and the properties of the phase space in dynami- 
cal systems of three degrees of freedom. One way to over- 
come this problem is to project the surface of the section to 
space with lower dimensions. In fact, here, we will apply 
the method introduced by Pfenniger (1984) (see also Revaz 
& Pfenniger 2001). We take sections in the plane y = 0, 
p y > of 3D orbits, whose initial conditions differ from the 
plane parent periodic orbits only by the z component. The set 
of the resulting four-dimensional points in (x, p x , z, p z ) phase 
space is projected on the (z, p z ) plane. If the projected points 
lie on a well-defined curve, we call it an 'invariant curve', 
then the motion is regular, while if not, the motion is chaotic. 
The projected points on the (z, p,) plane show nearly invari- 
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Figure 15: (a-b) Projections of the sections of 3D orbits with the plane y = when p y > 0. The set of the four-dimensional 
points (x, p x , z, p z ) is projected on the (z, p z ) plane. 



ant curves around the periodic points at z = 0, p z = 0, as 
long as the coupling is weak. When the coupling is stronger, 
the corresponding projections in the (z, p.) plane displays an 
increasing departure of the plane periodic point, up to mak- 
ing a direct orbit a retrograde one and vice versa. Here, we 
must define what one means by direct and retrograde 3D or- 
bit. If consequents in the (z, p.) section of the 3D orbit drop 
in one of the two domains of the corresponding section of 2D 
orbits at the same value of the Jacobi integral Ej, we can dis- 
tinguish between direct and retrograde motion. Orbits which 
visit both domains are intermittently direct or retrograde. 

Figures 15(a-b) show two typical (z,p z ) sections of 3D 
orbits, starting with initial conditions close to two different 
stable periodic points on the (x, p x ) phase planes of the 2D 
system. In order to obtain the results shown in Figure 15(a), 
we have taken the point (xo, Pm) = (0.25, 0) representing ap- 
proximately the position of the retrograde periodic point in 
the (x, p x ) phase plane of Figure 7(a). Similarly, in Figure 
15(b) we can observe the (z,p z ) projections near the retro- 
grade periodic point in the (x, p x ) phase plane of Figure 7(b). 
The position of the periodic point is (xo>P.<-o) = (0.25,0). 
Note that in both cases the numerical results indicate that for 
small values of zo the motion is regular, while for larger val- 
ues of zo the motion becomes chaotic. The integration time 
for each 'invariant curve' shown in Figures 15(a-b) is 3 x 
10 3 time units. Numerical calculations not given here sug- 
gest that the above method can be applied in all regular re- 
gions around each retrograde or direct stable periodic point 
of Figures 7(a-b), 8(a-b), 9(a-b), and 10(a-b). 

We must emphasize that the results presented in Figures 
15 (a-b) are rather qualitative and can be considered as an 
indication that the transition from regularity to chaos in 3D 
orbits occurs as the value of zo increases. In order to form 
a more complete and accurate view of the phase space in 
the 3D system, we computed a large number of 3D orbits 
(approximately 10 3 ) near each periodic point of the (x,p x ) 
phase planes of the 2D system, for different initial conditions 
(xo>P.io) and also for different values of zo- Our target was 
to determine the average minimum value of zo for which the 
nature of a 3D orbit changes from regular to chaotic. Table 
1 shows the value (zmm) near the direct and retrograde stable 
periodic points of Figures 7(a-b), 8(a-b), 9(a-b), and 10(a- 



b), for four different values of distance R between the centers 
of the two galaxies. From Table 1 we can induce three impor- 
tant results, (i) In every case (active or quiet primary galaxy) 
the minimum value of zo in the regions near the retrograde 
stable periodic points is always larger than the corresponding 
of the regions of the direct periodic points, (ii) If we compare 
each kind of galaxy regarding the nuclear activity, (active or 
quiet) for the same distance R, we observe that when the nu- 
cleus is not present at the galactic core (quiet - non-active 
primary galaxy), the 3D orbits can approach higher values 
of zo and remain regular. On the other hand, when the nu- 
cleus is present, the value of <z m j n } is smaller, (iii) As the 
distance between the centers of the galaxies increases, the 
minimum initial value of zo for which a 3D can remain reg- 
ular increases in both cases (active or quiet primary galaxy). 
This means that as the two galaxies are at large distances, 
their mutual interactions are weak enough and, therefore, the 
majority of the 3D orbits are ordered. Moreover, at large 
distances, one can conclude that the main factor responsible 
for the observed chaotic motion is the nuclear activity of the 
nucleus in the core of the primary galaxy. 

So far, we have seen that 3D orbits with initial condi- 
tions (xo,p^o,Zo)> sucn that (x ,p. r o) is a point in the chaotic 
regions of the 2D system, for all permissible values of zo are 
chaotic. On the other hand, the nature (ordered or chaotic) 
of 3D orbits with initial conditions (x ,p^o,Zo), such that 
(->co,P.io) is a point in the regular regions around the stable 
direct or retrograde periodic points of the 2D system, de- 
pends on the particular value of zo> as we can see in Table 
1. We did not feel that it was necessary to try to define the 
values of <z m i n > for each regular region of the 2D system cor- 
responding to secondary resonances which are represented 
by sets of multiple small islands of invariant curves in the 
(x, p x ) phase planes shown in Figures 7(a~b), 8(a~b), 9(a-b), 
and 10(a-b). Numerical results indicate that 3D orbits with 
initial conditions (xo,p A o,Zo), such that (x Q ,p M ) is a point in 
the regular regions corresponding to secondary resonances of 
the 2D system, remain regular for (z m i„) 0.052, while for 
larger values of z , they change their character from regular 
to chaotic. 

In Figures 16(a-h) we present eight typical orbits of the 
3D dynamical system. Figure 16(a) shows a regular quasi- 
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Figure 16: (a-h) Eight representative regular orbits of the 3D dynamical system. The initial conditions and more details 
regarding the values of all the other parameters are given in the text. 



periodic orbit circulating around the primary galaxy, with ini- 0, while the value of p y0 is always found from the Jacobi in- 
tial conditions x = 0.167, y = 0, zo = 0.01, and p M = p- = tegral (10). In Figure 16(b), we observe a 3D periodic orbit 
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Table 1: Average Value of the Minimum zo Near the 
Direct and Retrograde Periodic Points, for Four Differ- 
ent Values of the Distance R 



Distance 


Case 


Region 




R - 1.5 


rVCLlVC 


Direct 


- 


Retrograde 


0.0792 




Direct 


- 


Retrograde 


0.0802 


R - 2.0 


A r* 1 1 \ 7£» 
r\CLlVC 


Direct 


0.0865 


Retrograde 


0.0912 


Quiet 


Direct 


0.0896 


Retrograde 


0.0964 


R - 2.5 


Active 


Direct 


0.0983 


Retrograde 


0.1025 


Quiet 


Direct 


0.1012 


Retrograde 


0.1107 


R - 3.0 


Active 


Direct 


0.1098 


Retrograde 


0.1157 


Quiet 


Direct 


0.1185 


Retrograde 


0.1214 



characteristic of the 1 :2 resonance. This orbit has initial con- 
ditions x = -0.1895, y = 0, zo = 0.02, and p t0 = p z0 = 0. 
Figure 16(c) shows a 3D periodic orbit with initial conditions 
Xq = 0.5788, y = 0, zo = 0.01, and p M = p z0 = 0, circu- 
lating around the center of the primary galaxy. In Figure 
16(d), we see a quasi-periodic orbit, with initial conditions 
Xq = 0.4295, y = 0, z = 0.015, and p M = p z0 = 0. 
A 3D resonant orbit with initial conditions Xo = 0.7445, 
Vo = 0, zo = 0.01, and p M = p zQ = 0, moving around the 
galactic center, is shown in Figure 16(e). In Figure 16(f) we 
observe a 3D periodic orbit circulating around the primary 
galaxy. This orbit has initial conditions x = -0.9854, y = 
0, zo = 0.01, and p^ = p z o = 0. In Figure 16(g) we observe 
a complicated 3D resonant periodic orbit with initial condi- 
tions x = 0.673, y = 0, zo = 0.005, and p v0 = p z o = 0. 
Figure 16(h) shows a 3D chaotic orbit with initial conditions 
Xo = -0.31, jo = 0, zo = 0.1, and p M = p z0 = 0. This orbit 
goes arbitrarily close to the primary galaxy and it is deflected 
to higher values of z, on approaching the dense and massive 
nucleus. We must note that in all 3D orbits shown in Figures 
16(a-h), the initial conditions (xo,Pm) and the values of all 
the other parameters (M d , M„, R) are as in the corresponding 
2D orbits shown in Figures 13(a-h). Moreover, we observe 
that all regular 3D orbits shown in Figures 16(a-g) stay near 
the galactic plane and, therefore, support the disk structure 
of the primary galaxy. The numerical integration time for all 
3D orbits shown in Figures 16(a-h) is 200 time units. 

5 A Theoretical Approach 

In this section, we shall present some theoretical arguments, 
together with elementary numerical calculations, in order to 
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Figure 17: A plot of the total force as a function of the 
distance R, when the primary galaxy is active or quiet. 



explain the numerically obtained relationships given in Fig- 
ures 1 1, 12, and 14. Moreover, we will quote different kinds 
of theoretical techniques in an attempt to explain and study 
deeper the dynamical structure of the 2D or 3D Hamiltonian 
system. 

Since the potential of the nucleus is integrable with spher- 
ical symmetry, the observed chaotic phenomena in our dy- 
namical system should derive mainly from the total \F,\ force. 
Figure 17 shows a plot of the \F,\ force as a function of the 
distance R between the two galaxies at the point P c (xo, yo, zo) = 
(-0.12,0.02,0.01). This point is very close to the center of 
the primary galaxy, where the chaotic phenomena are more 
prominent. One may observe, from the plot shown in Fig- 
ure 20 that the two curves, which correspond to two differ- 
ent cases (active and quiet primary galaxy), have almost the 
same pattern. In both cases, the |F,| force decreases, tend- 
ing asymptotically to zero, as the value of the distance R 
increases. Here, we must point out that the values of the 
|F,| force are larger when the primary galaxy is active and, 
therefore, the corresponding curve is above the one corre- 
sponding to the quiet galaxy. Thus, we conclude that when 
a massive and dense nucleus is present in the center of the 
primary galaxy, the extent of the chaotic orbits is larger than 
in the case in which the primary galaxy is quiet. One can 
see that the pattern between the force and the distance R 
shown in Figure 17 and obtained theoretically is very sim- 
ilar to those connecting the chaotic percentage A% or the 
LCE with the distance R shown in Figures 11, 12, and 14, 
respectively, which have been obtained numerically. Thus, 
the numerically obtained results given in Figures 11, 12, and 
14 regarding the evolution of the chaos with respect to the 
distance R can be justified and explained, in a way, theoret- 
ically through the relation between the strength of the total 
force and the distance R depicted in Figure 17. 

Furthermore, it would be of particular interest to study 
the structure of the velocity profile that is the plot of the 
total velocity of a test particle (star), v = ^jx 2 + y 2 + z 2 as 
a function of time, for ordered and chaotic motion. Figure 
18(a) shows the velocity profile for a time period of 500 time 
units corresponding to the regular 3D orbit shown in Figure 
16(a). Here, the velocity profile is quasi-periodic and the 
maximum value of the velocity is about 450 kms~'. This 
suggests that the regular 3D motion occurs in small veloci- 
ties. Figure 18(b) shows the velocity profile for a time pe- 
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Figure 19: A plot of the L- component of the total angular momentum versus time for (a) the regular 3D orbit shown in Figure 
16(a) and (b) the chaotic 3D orbit shown in Figure 16(h). 



riod of 500 time units corresponding to the chaotic 3D orbit 
shown in Figure 16(h). In this case, we can discuss two as- 
pects. The first is that the velocity obtains high values up to 
790 km s -1 and the second is that the velocity profile appears 
to be highly asymmetric, displaying abrupt changes and large 
deviations between the maxima and also between the min- 
ima in the [v - t] plot. Therefore, we may conclude that the 
chaotic 3D motion occurs in high and abruptly changing ve- 
locities. 

The physical parameter playing an important role in the 
orbital behavior of the stars is the L. component of the total 
angular momentum. From our previous experience, we know 
that on approaching a dense and massive nucleus, low angu- 
lar momentum stars are scattered off the galactic plane, dis- 
playing chaotic motion (Caranicolas & Innanen 1991; Caran- 
icolas & Papadopoulos 2003; Caranicolas & Zotos 2010; Zo- 
tos 201 la, 201 lc). Of course, here in 3D phase space, things 
are more complicated than in axially symmetric dynamical 
models, where the L- component is conserved. As the mo- 
tion takes place in a rotating non-axially symmetric system, 
the L- component is not conserved and it is given by 

L z = xy-xy-Q p [x 2 + y 2 ) . (15) 
Nevertheless, we can compute numerically its mean value 



(Lj> using the formula 

1 " 

<L ; >=-VL-„ (16) 
n to 

Our numerical calculations suggest that the chaotic orbits 
have low values of (L-), while regular orbits obtain high val- 
ues of (L z ). Figure 19(a) shows a plot of the evolution of 
the L z component with the time for the regular orbit of Fig- 
ure 16(a). In this case, we observe that L- is nearly a quasi- 
periodic function of time, while (L ; ) = -0.12271. Figure 
19(b) is similar to Figure 19(a) but for the chaotic orbit shown 
in Figure 16(h). Here one can see abrupt changes in L, dur- 
ing the chaotic motion, while for this chaotic orbit we have 
(L z ) = -1.24662. In both cases, the time interval of the nu- 
merical integration is 500 time units, while n = 10 4 . 

In what follows we shall present a semi-theoretical anal- 
ysis, in order to give a more detailed picture of the structure 
of the dynamical system and its behavior. The forces act- 
ing on a test particle along the x and y axes are given by the 
equations 

M n (x-xj) M A (x-Xj) 
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Figure 21: (a) A 3D plot of the value of the F, force on the (x,z) plane and (b) contours of the projections F- = const on the 
{x,z) plane. Lighter colors indicate higher values of F-. For positive values of z the F z force is negative, while for negative 
values of z the F z force is positive. 
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Figure 20: Contours of the F x = const (elliptic shaped) 
together with the contours F y = const (figure-eight-shaped). 
Details are given in the text. 

M s (x - x 2 ) , 

- : + 0. x - 2Cl p y, 

p M n y M&y 

~{t l+ cf 2 ~[bi + l j ll+ (a + h)f 2 

- ^f+nly + 2n p x. (17) 

It is obvious from Equations (17) that the strength of both 
forces increases as the mass of the nucleus or the disk in- 
creases or their scale lengths decrease. Figure 20 shows the 
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Figure 22: A plot of the total velocity v(x) as a function 
of the distance x, for the two cases (active and quiet primary 
galaxy). 



contours of F x = const together with the contours F y = 
const. The values of all other parameters are as in Figure 
2. The contours of F x = const look like ellipses, while the 
contours of F y = const have a figure-eight shape. Inside the 
first 'ellipse' we have F x > 0, while inside the second 'el- 
lipse' we have F x < 0. On the other hand, inside the lower 
figure-eight curve we have F y > 0, while inside the upper 
figure-eight curve we have F y < 0. Looking carefully near 
the galactic center, we observe that there are areas where F x 
is positive and F y is negative at the same time and vice versa. 
There are also areas where both forces are positive or nega- 
tive at the same time. Therefore, we conclude that near the 
center of the primary galaxy, there are strong attractive and 
repulsive forces acting on the test particle (star). These forces 
are responsible for the chaotic scattering of the star near the 
galactic center, leading to chaotic motion. 

At this point, we can expand the semi-theoretical anal- 
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ysis described above, in order to give a more detailed and 
complete picture of the structure of the 3D dynamical sys- 
tem and its behavior. The force acting on a test particle along 
the z axis is given by the equation 



M d (a + r z )z 



3/2 



(18) 



where r z = yh 2 + z ■ It is obvious from Equation (18) that 
the strength of the F, force increases as the mass of the nu- 
cleus or the disk increases or their scale lengths decrease. 
Figure 21(a) shows a 3D plot of the value of the F z force on 
the (x, z) plane. Figure 21(b) depicts the contours of the pro- 
jections F z = const on the (x, z) plane. We observe that for 
positive values of z, the F z force is negative, while for nega- 
tive values of z, the value of the F, is positive, near the center 
of the primary galaxy. Therefore, lighter colors on the lower 
half part of the (x, z) plane indicate higher values of the F z 
force, while darker colors on the upper half part of the same 
plane indicate lower values of the F z force. One can observe 
that near the active galactic core which hosts a dense and 
massive nucleus, the test particle experiences a very strong 
F z force. The values of all the parameters in Figures 21(a-b) 
are as in Figure 2. 

Last but not least, we also investigate the behavior of 
the total velocity near the galactic core as a function of the 
distance x. To do that, we consider the limiting curve, that 
is, the curve containing all the invariant curves on the (x, p x ) 
phase plane. This can be obtained if we set y = z = p y = 
p z = in Equation (10), yielding 



1 



l 



(19) 



where we have set p x = v because at the limiting curve the 
p x velocity is the total velocity (see also Papadopoulos & 
Caranicolas 2006). In Figure 22, we observe two plots of 
the total velocity v as a function of the distance x, derived 
using relation (19). The values of all the parameters are as 
in Figure 2. The dashed line corresponds to the case of the 
quiet primary galaxy, while the solid line corresponds to the 
case when the primary galaxy hosts a nucleus in its core and, 
therefore, is active. We see that the velocity is about the same 
for all values of x, except near the galactic center, that is 
when x =s -0.13, where higher velocities correspond to the 
active galactic center, hosting massive and dense nucleus. 

One of the factors responsible for the chaotic motion and 
other resonance phenomena, such as islandic motion corre- 
sponding to sets of multiple islands of invariant curves or 
invariant tori, are the several inner Lindblad resonances 



ci„ = ci- 



-K, 



(20) 



where CI and k indicate the circular and the epicycle fre- 
quency of the star, respectively, while m and n are integers. 
The main resonances for the two cases (active and quiet pri- 
mary galaxy) together with the corresponding resonance radii 
n and r 2 , when R = 2.35 and Cl p = 0.419146, are given in 
Table 2. 

Figure 23(a) shows a plot of the curves fi - iik/iti in the 
case where the primary galaxy is active, as a function of the 
radius r, when R = 2.35. We choose this particular value 
of the distance between the two galaxies because at this dis- 
tance we can observe a large variety of Lindblad resonances. 



Table 2: Radii of Lindblad Resonances When/? = 2.35 
and Q p = 0.419146, for Both Active and Quiet Pri- 
mary Galaxies 



Case 


Resonance 


Region 


r\ 




Active 


lip — 11 — LKj J 


Direct 


0.0149 


- 


Retrograde 


0.0142 


0.2385 


lip — 11 — 


Direct 


0.0246 


- 


Retrograde 


0.0217 


0.3548 


CI - CI - 5*r/8 

lip — i£ .Jft/O 


Direct 


0.0127 


0.3886 


Retrograde 


0.0122 




Clp = Q - 7k/9 


Direct 


0.0257 




Retrograde 


0.0238 


0.3912 


Quiet 


Cl p =Cl- 2k/5 


Direct 


0.0123 


0.5527 


Retrograde 


0.0114 


0.1342 


Clp = CI- 3k/4 


Direct 


0.0120 




Retrograde 


0.0119 


0.2351 


Cl p = £1- 3k/7 


Direct 


0.0095 


0.6782 


Retrograde 


0.0091 




Clp = CI- 4k/9 


Direct 


0.0207 




Retrograde 


0.0196 


0.2608 



The numbers 1, 2, 3, and 4 indicate the curves CI - 2k/3, 
CI — 3k/4, Q. - 5/c/8, and CI - Ik/9, respectively. The straight 



lines are the curves Cl„ 



t0.419146. The values of all the 



other parameters are as in Figure 2. Figure 23(b) is similar to 
Figure 23(a) but when the primary galaxy is quiet. Here, the 
numbers 1, 2, 3, and 4 indicate the curves C1 — 2k/5, CI- 3k/ 4, 
C1-3k/7, and C1-4k/9, respectively. As we can see, there are 
a considerable number of resonance radii for both the direct 
and retrograde orbits in both cases. Details regarding the res- 
onances and the resonance radii can be obtained from Table 
2. In other words, all the Lindblad resonances given in Table 
2 are also responsible for the chaotic motion in the primary 
galaxy. It is also interesting to note that the above resonances 
produce large chaotic regions for small values of the distance 
R, while for larger values of R (see Figures 9a-b) the chaotic 
regions are extremely small, although the resonance radii are 
still present. This means that in this case the distance be- 
tween the two galaxies precedes the Lindblad resonances. 

In all cases, the primary galaxy is assumed to have a fixed 
value of total mass 



M P = M d + M n = 2.08, 



(21) 



where < M n < 0.08. For a constant distance R, the F x , F y , 
and F, forces, per unit mass, at an arbitrary but fixed point 
Po(x ,y ,zo) near the center of the primary galaxy are linear 
functions of the mass of the nucleus M n , of the form 



F x = 




-kiM D , 




Fy = 




- k 2 M n , 




Fz = 


-c 3 


- k 3 M n , 


(22) 



where (ci , c 2 , c 3 ) and (k t ,k 2 ,k 3 ) are positive constants, while 
relation (21) has also been taken into account. Figure 24(a) 
shows the total \F,\ force as a function of M n , at P (x ,y Q , zo) = 
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Fi gure 23; The curves Q — nK/m versus r, when R = 2.35 and Q p = 0.419146, when (a) the primary galaxy is active and (b) 
the primary galaxy is quiet. 




0.02 0.04 0.06 0.08 0.02 0.04 . OS 0.08 




Figure 24: (a-d) Evolution of different dynamical quantities of the system as a function of the mass of the nucleus M n . Details 
are given in the text. 



(-0.12,0.02,0.01), when R = 1.5. We observe that the to- 
tal force increases rapidly as the mass of the nucleus M n in- 
creases. The presence of the nucleus increases the velocity 
near the center of the primary galaxy. This can be easily 
shown if we consider the limiting curve, that is, the curve 
enclosing all the invariant curves in the (x,p x ) phase plane. 
This is obtained if we set y = in the Jacobi integral (13). 
Thus, we have 



The curve 



n 2 = 2[£y-(D,(x)]-p2>o. 



/(*, p,) = 2 [£ y -*,(*)] -/*[ = 



(23) 



(24) 



is a curve on the (x, p.,) phase plane, called the limiting curve 
or the zero velocity curve (ZVC). The velocity u. v = p x at the 
limiting curve is the total velocity v. At a fixed point x = x , 
the total velocity is 



2Ej + const + c 4 (xo - Xi) + 



1/2 



(25) 



where c 4 < 1 is a constant, while 



yjb 2 + (jco - xi) 2 + (a + h) 2 , 
Jixo-xrf + c*. (26) 
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As A\ > A 2 , the maximum value of the velocity is met when 
x Q = x, (see also Figures 7a-b). It is evident from Equations 
(25) and (26) that v increases as M„ increases or c„ decreases. 
This indicates that high velocities are expected in the central 
regions of galaxies with dense and massive nuclei. Figure 
24(b) shows a plot of v as a function of M n when R = 1 .5 . We 
see that the velocity increases as the mass of the nucleus M n 
increases. In Figure 24(c) we observe a plot of the rotational 
velocity © of the primary galaxy, as a function of the mass of 
the nucleus M n , when r = r = 0.01. The exact function is 



&(M n ) = r . 



M n 



\(rl+clY {V + rl + {a + h?) 



3/2 ' 



(27) 



From Figure 24(c) we can derive two basic conclusions, (i) 
In galaxies with massive and dense central concentrations, 
the rotational velocity has and retains high values near the 
galactic center, (ii) The more massive and dense the nucleus, 
the higher is the value of the rotational velocity 0. Figure 
24(d) shows the evolution of the mass density p of the pri- 
mary galaxy, as a function of the mass of the nucleus M n 
at the point P (x ,y ,zo) = (-0.12,0.02,0.01). It is evident 
that as the galaxy evolves and a nucleus forms in its center, 
we have a constant and rapid increase of the mass density. 
The relationships given in Figures 24(a-d) are of great im- 
portance and will prove very useful in the following section, 
where we will study the evolution of our galactic system us- 
ing a time-dependent model. 



6 Evolution of the Orbits in the 3D 
Time-Dependent Model 

Let us now follow the evolution of 3D orbits, as mass is trans- 
ported from the disk of the quiet galaxy to its center. By this 
procedure, a massive and dense nucleus is developed in the 
central regions of the primary galaxy. The mass transport is 
linear following the set of equations: 



M nf (/) = M m + kt, 
M if (t) = M di - kt, 



(28) 



where M ni = and M di = 2.08 are the initial values of the 
mass of the nucleus and the disk respectively, while k is a 
positive parameter. We assume that the linear rate described 
by relations (28) is slow compared to the orbital period of the 
binary system and therefore it is adiabatic. This is true be- 
cause the mass transportation period is 100 time units, while 
the orbital period is about three orders of magnitude smaller. 
It is also assumed that the transportation stops when the mass 
of the nucleus in the galactic core takes the value M n f = 0.08. 
It is well known that the shape of the 3D orbits sometimes is 
inconclusive or even misleading. In order to overcome this 
drawback, we have decided to use a more accurate method 
such as the LCE. The main advantage of this method is that 
it uses certain and objective numerical thresholds, beyond 
which we can distinguish between ordered and chaotic mo- 
tion. On the other hand, this method is very time-consuming 
as it needs time intervals of numerical integration of the order 
of 10 5 time units, in order to provide reliable and definitive 
results regarding the nature of a 3D orbit. But this is a 'price' 
we can afford to pay. Additionally, we shall use a fast but 



qualitative method, that is the profile's pattern of the total ve- 
locity, so as to characterize a 3D orbit. In the following, we 
present the evolution of four different 3D orbits, as the to- 
tal mass distribution of the dynamical system changes with 
time, following the set of equations (28). For all orbits shown 
in Figures 25-28, the initial value of the Hamiltonian (10), at 
t = 0, is Eji = -2.6, while k = 0.0008. 

In Figure 25(a) we can see the evolution of the LCE for 
a 3D orbit and for a time period of 10 5 time units, as the 
galaxy evolves following the set of equations (28). The ini- 
tial conditions are x = 0.12, y = 0, zo = 0.02, p M = 0, and 
p z0 = 0, while the value of p y0 is found from the Hamiltonian 
(10) in all cases. The values of all the other parameters are 
as in Figure 7(b). When t = 100 time units, the mass of the 
developed nucleus in the core of the galaxy is M nf = 0.08 
and the evolution stops. The value of the Hamiltonian is now 
Ejf = -2.60012. The profile of the LCE clearly indicates 
that this orbit starts as a regular and remains regular during 
the galactic evolution. Figure 25(b) shows the velocity pro- 
file for the same 3D orbit for a time period of 200 time units. 
The pattern is almost quasi-periodic and therefore the orbit 
remains regular during the galactic evolution. 

On the other hand, things are quite different in Figure 
26(a). This orbit has initial conditions x = -0.44, y = 0, 
Zo = 0.01, p x0 = 0, and p z0 = 0. The values of all the 
other parameters are as in Figure 7(b). When t = 100 time 
units, the mass of the developed nucleus in the core of the 
primary galaxy is M nf = 0.08 and the mass transportation 
stops. The value of the Hamiltonian is now Ejf = -2.60057. 
In this case, the profile of the LCE clearly indicates that this 
orbit starts as a chaotic 3D orbit and remains chaotic dur- 
ing the galactic evolution. Figure 26(b) shows the velocity 
profile for the same 3D orbit, for a time period of 200 time 
units. The pattern is highly asymmetric, with a large number 
of abrupt peaks, large deviations between the maxima and 
also significant deviations between the minima. Therefore, 
we conclude that this orbit remains chaotic during the mass 
transportation. 

In Figure 27(a) we observe the evolution of the LCE 
of a 3D orbit with initial conditions x = -0.196, yo = 0, 
Zo = 0.01, p x0 = 0, and p z0 = 0. The values of all the other 
parameters are as in Figure 7(b). After a time interval of 100 
time units the mass of the nucleus in the core of the galaxy 
is M nf = 0.08 and the galactic evolution stops. The Hamil- 
tonian settled to the value Ejf = -2.60034. The profile of 
the LCE given in Figure 27(a) shows that the orbit starts as 
a regular 3D orbit, but after the galactic evolution it becomes 
chaotic. It is evident that if mass transport were not present, 
the orbit would have remained regular. The presence of the 
nucleus in the core of the primary galaxy has changed the 
character of the 3D orbit from regular to chaotic. The same 
result can be obtained by looking at the profile of the veloc- 
ity of the 3D orbit, as shown in Figure 27(b). For the first 
100 time units the pattern is quasi-periodic, indicating reg- 
ular motion, but after that time interval it becomes highly 
asymmetric, leading to the conclusion that the orbit has be- 
come chaotic. 

Figure 28(a) depicts the evolution of the LCE of a 3D 
orbit with initial conditions x = 0.74, yo = 0, zo = 0.01, 
Pxo = 0, and p z0 = 0. The values of all the other parame- 
ters are as in Figure 7(b). After a time interval of 100 time 
units, the mass of the nucleus in the core of the galaxy is 
M nf = 0.08 and the mass transportation stops. The Hamilto- 




Figure 25: (a) Evolution of the LCE of the 3D orbit in the time-dependent model, following relations (28), and (b) the 
corresponding velocity profile. The orbit starts as a regular and remains regular during the mass transportation. 

W (b) 




Figure 26: (a) Evolution of the LCE of the 3D orbit in the time-dependent model, following relations (28), and (b) the 
corresponding velocity profile. The orbit starts as a chaotic and remains chaotic during the mass transportation. 

(a) (b) 




Figure 27: (a) Evolution of the LCE of the 3D orbit in the time-dependent model, following relations (28), and (b) the 
corresponding velocity profile. The orbit starts as a regular but after 100 time units; when the galactic evolution stops it becomes 
chaotic. 



nian settled to the value E Jf = -2.60021. The profile of the chaotic 3D orbit, but after the galactic evolution it becomes 
LCE given in Figure 28(a) indicates that this orbit starts as a an ordered one. Therefore, it is evident that if mass trans- 
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Figure 28: (a) Evolution of the LCE of the 3D orbit in the time-dependent model, following relations (28), and (b) the 
corresponding velocity profile. The orbit starts as a chaotic but after the galactic evolution, it changes its nature to regular. 



portation were not present, the orbit this time would have 
remained chaotic. In this case, the presence of the nucleus in 
the core of the galaxy has changed the nature of the 3D or- 
bit from chaotic to regular. A similar result can be obtained 
from the profile of the velocity of the 3D orbit, as shown in 
Figure 28(b). For the first 100 time units the pattern is highly 
asymmetric, indicating chaotic motion, but for the rest 100 
time units, it becomes almost quasi-periodic, leading to the 
conclusion that the orbit has now become regular. 

The above analysis reveals that as the galaxy tends to 
create the central nucleus, the character of most 3D regular 
orbits, starting near the center of the primary galaxy, changes 
their nature from regular to chaotic. At the same time, the ve- 
locity shows significant changes and increases as the mass of 
the nucleus increases. In other words, as the quiet galaxy 
shows this active face, the nature of motion presents ma- 
jor alterations. Using the above procedure, we have tested 
a large number of 3D orbits (approximately 10 3 ) in the time- 
dependent model, describing the formation of a massive and 
dense nucleus in the core of the primary galaxy, when mass 
is transported from the disk. Our numerical results indi- 
cate that the character of the 3D orbits can change either 
from regular to chaotic and vice versa or not change at all, 
as the mass is transported in order to create the nucleus in 
the central region of the disk galaxy. In particular, from 
the sample of the 10 3 tested orbits in the case of the time- 
dependent model, we conclude that 61% of the orbits altered 
their nature from regular to chaotic, 22% remain chaotic, 
15% remain regular, and only 2% changed their character 
from chaotic to regular. Thus, it is evident that the forma- 
tion of a massive nucleus leads the majority of the orbits of 
the 3D system to become chaotic. Whether the nature of 
a 3D orbit will change or not during the galactic evolution 
described by the set of equations (28) depends strongly on 
the initial conditions (xo,p x o,zo) of each orbit. Moreover, 
as the change of the value of the Hamiltonian (10) is negli- 
gible (A£j = \Ejf - Ejj\ a 10~ 4 ), we can say that the phase 
space is transformed to itself during the galactic evolution 
and, therefore, the orbits can be considered isoenergetic. In 
order to make this statement more clear, we present an exam- 
ple. If we suppose that the evolving time-dependent model 
was describing the 2D system, then as the change of the 



Hamiltonian is negligible, we could say that in the phase 
plane of the quiet galaxy, shown in Figure 7(b), chaotic re- 
gions would appear in the central area and it would be trans- 
formed to the phase plane shown in Figure 7(a), as the nu- 
cleus is formed in the core of the primary disk galaxy. 

7 Discussion and Conclusions 

In the present article, we have constructed a 3D gravitational 
model, in order to study the character of motion in a binary 
galactic system. In particular, we present a dynamical model 
which is composed of a primary disk galaxy and a small 
satellite companion. The motion of a test particle (star) in 
the gravitational field of this galactic system was studied and 
various techniques were used to identify regular and chaotic 
motion. The galaxies were assumed to be coplanar and or- 
biting each other in the same plane on a circular orbit. In 
particular, this research is indeed an extension to three di- 
mensions of the work of Caranicolas & Innanen (2009), who 
studied a similar two-dimensional binary system of interact- 
ing galaxies. 

A binary system of interacting galaxies is very complex 
and, therefore, we need to assume some necessary simplifi- 
cations and assumptions, in order to be able to study the or- 
bital behavior of such a complicated stellar system. Thus, our 
model is simple and contrived, in order to give us the ability 
to study different aspects of the dynamical model. Never- 
theless, contrived models can provide an insight into more 
realistic stellar systems, which unfortunately are very diffi- 
cult to be studied if we take into account all the astrophysical 
aspects. Here, we must point out the main restrictions and 
limitations of our gravitational model: (1) The two galaxies 
(the primary and the satellite) are assumed to be coplanar, or- 
biting each other in the same plane on circular orbits. (2) Our 
dynamical model only deals with the non-dissipative com- 
ponents of the galaxies, stars, or possibly dark matter parti- 
cles. (3) The potentials we use are rigid and do not respond 
to the evolving density distribution in a more realistic way. 
This is because our gravitational model that describes the bi- 
nary system of the interacting galaxies is not self-consistent. 
Thus, the two galactic centers remain stationary in the rotat- 
ing frame chosen. Self-consistent models are usually used 
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when conducting A'-body simulations. Obviously, this is out 
of the scope of the present paper. Once again, note that the 
above restrictions and limitations of our model are necessary; 
otherwise it would be extremely difficult, or even impossi- 
ble, to apply the extensive and detailed dynamical study pre- 
sented in this study. A similar gravitational model with the 
same limitations and assumptions was used in Zotos (2012b) 
in order to study the motion in a binary system of interacting 
galaxies, where the second galaxy is no longer treated as a 
point-mass satellite companion. Thus, we may conclude that 
the same setup has been successfully applied to describe a 
much more complex binary galactic system. 

In our study, two different cases were investigated: the 
time-independent model and the time-evolving model, that 
is the case when mass is transported adiabatically from the 
disk of the galaxy to its center, forming a massive and dense 
nucleus. Our numerical calculation indicates that there are 
several factors responsible for the observed chaotic motion 
in the time-independent model: (i) the galactic interaction, 
(ii) the galactic activity, that is the presence of the nucleus, 
and (iii) the Lindblad resonances. Furthermore, the presence 
of the nucleus increases the velocities near the central region 
of the primary galaxy. The value of the velocity depends 
on the mass of the nucleus and also on the value of its scale 
length. Regular motion corresponds to low central velocities, 
while chaotic motion is characterized by high velocities. All 
the above observations strongly indicate that in the centers 
of active galaxies, chaotic motion in high velocities is ex- 
pected. On the other hand, it was found that the two interact- 
ing galaxies, for large values of the distance R between their 
centers, do not present chaotic motion, when the nucleus is 
not present in the core of the primary galaxy. The results 
of this work regarding the nature of the orbits and also the 
factors that affect it or change it are very similar to the corre- 
sponding outcomes obtained in Zotos (2012b) and, therefore, 
verify in a way our previous conclusions. 

We have started our investigation from the Hamiltonian 
system of two degrees of freedom (2D). Our numerical cal- 
culations indicate that in this case, a large part of the phase 
plane is covered by chaotic orbits, while the regular regions 
are confined mainly near the stable retrograde periodic point. 
The chaotic area is larger when the primary galaxy possesses 
a dense and massive nucleus or when the distance R between 
the two galaxies is small. It was also found that the velocity 
near the central region of the galaxy increases significantly 
when the nucleus is present. For lower values of the total 
energy and small values of the distance R, no chaos is ob- 
served when the nucleus is absent, while a small chaotic re- 
gion appears in the presence of the nucleus. This means that 
low-energy stars are in chaotic orbits near the centers of ac- 
tive galaxies with a suitable satellite companion in a circu- 
lar orbit, while in quiet galaxies they are not. Furthermore, 
the velocity near the center rises to high values, even for 
smaller values of the energy. The above discussion strongly 
indicates that in the centers of active galaxies, chaotic mo- 
tion and high velocities are expected. This fact, combined 
with outcomes from previous works, shows that the major- 
ity of orbits in galaxies with dense and massive nuclei are in 
chaotic orbits (see Caranicolas & Innanen 1991; Caranicolas 
& Papadopoulos 2003; Zotos 2012a). This seems reasonable 
because theoretical results show that the nuclear force near 
the center increases linearly as M n increases. Furthermore, 
the velocity near the galactic center strongly depends on the 



scale size of the nucleus c n . The smaller the c n , the higher is 
the velocity near the central region of an active galaxy. 

To explore and understand the nature of orbits in the 3D 
dynamical system, we have used our knowledge obtained 
from the study of the 2D system. Of particular interest was 
the determination of the regions of initial conditions in the 
{x,p x , z), p y > 0, (y = p z = 0) phase space that produce reg- 
ular or chaotic 3D orbits. As the value of p, was found from 
the Jacobi integral (10), we have used the same value of Ej 
as in the 2D system and took initial conditions {xQ,p M ,zo) 
such that (x ,p. v o) lies in the chaotic regions of the 2D sys- 
tem. It was found that the motion is chaotic for all permis- 
sible values of zo- On the other hand, when (x Q ,p xQ ) was 
inside a regular region around the direct and retrograde peri- 
odic points, the corresponding 3D orbits are regular for small 
values of zo, while for larger values of zo the orbits become 
chaotic. The particular values of zo for which the transition 
from regularity to chaos in 3D orbits occurred were different 
for each regular region of the 2D system. Of particular in- 
terest are the results given in Table 1, where we define the 
average minimum value of zo near the direct and retrograde 
periodic points of the primary galaxy, where the transition 
from regularity to chaos occurs. 

In this research, we have used different kinds of theoret- 
ical techniques in an attempt to explain and justify the dy- 
namical structure of the 2D or 3D Hamiltonian system. An 
important role is played by the L- component of the test par- 
ticle's angular momentum. It was found that the values of 
(L-) for regular 3D orbits are larger than those for chaotic 3D 
orbits. Thus, the L z component of the angular momentum is 
a significant dynamical parameter connected with the regular 
or chaotic character of orbits in both 2D and 3D dynamical 
systems. In order to estimate the degree of chaos in the 2D as 
well as in the 3D dynamical system, we have computed the 
average value of the LCE for a large number of orbits with 
different initial conditions in the chaotic regions in each case, 
for a time period of 10 5 time units. Our numerical results in- 
dicate that the degree of chaos in the 3D binary system of 
interacting galaxies is smaller than that in the 2D system. 

It was of great interest to follow the evolution of orbits 
as nucleus formed in the central region of the primary disk 
galaxies. In this procedure, mass was transported from the 
disk to the nucleus and therefore the quiet galaxy became 
gradually active, following Equation (28). We observed that 
the final character of the 3D orbits strongly depends on the 
particular initial conditions (xo, p x o,Zo)- Thus, regular orbits 
can change to chaotic and vice versa, or maintain their char- 
acter (ordered or chaotic) during the nuclear formation and 
since the nucleus has been formed. It was observed that 
a number of regular quasi-periodic orbits, starting near the 
central region of the primary galaxy, become chaotic. This 
can be seen in a total velocity versus time plot, where the 
asymmetric profile of the total velocity is in agreement with 
observational data (see Grosb0l 2002), where an increase of 
the stellar velocity is expected, in regions with significant 
chaoticity. Moreover, observations show that an asymmet- 
ric velocity profile indicates chaotic motion. An interesting 
question is whether interactions are essential in order to trig- 
ger the mass transportation and therefore the galactic activity. 
This question remains unsolved as astronomers have found 
binary systems with disk galaxies (see Letawe et al. 2006) 
showing no signs of interaction at all but harboring active 
nuclei. 
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Fi gure 29 1 A real image of the binary galactic system com- 
posed of the primary galaxy NGC 5829 and its small satellite 
companion IC 4526. 



Our theoretical dynamical analysis of the binary galactic 
system can now be compared with observational data derived 
from the M51-type binary galactic systems. These systems 
consist of a large, primary disk galaxy and a smaller satel- 
lite companion galaxy moving in a circular orbit around their 
common mass center (see Klimanov & Reshetnikov 2001). 
We apply our results to the binary system composed of the 
primary galaxy NGC 5829 and its small satellite IC 4526. 
Figure 29 depicts a real image of this binary system. The 
primary galaxy has an absolute magnitude of M v = -20.66. 
In order to estimate the disk's radius, we shall use the empir- 
ical relation 

M v = -6 log A + 7.14, (29) 

which connects the absolute magnitude to the projected ma- 
jor axis of the primary galaxy A in pc (see Bowers & Deem- 
ing 1984). From relation (29), we obtain the value A = 43 
kpc. If we adopt the values of the parameters given in Figures 
7(a-b), then one should expect a large number of the orbits 
to be chaotic in this binary system of interacting galaxies. 
We can also observe that the orbit of the satellite is close to 
the boundary of the galactic disk. This result is in agreement 
with the observational data as well. 

When studying galactic dynamics, it is invariably found 
that the stellar rotation velocity remains constant, or flat, with 
increasing distance away from the galactic center (see Can- 
tinella, Giovanelli, & Haynes 2006). This result is highly 
counterintuitive since, based on Newton's law of gravity, the 
rotational velocity would steadily decrease for stars farther 
away from the galactic center. By this particular argument, 
the flat rotational curves seem to imply that each galaxy must 
be surrounded by significant amounts of dark matter. It has 
been postulated and generally accepted that the dark matter 
would have to be located in a massive halo enshrouding each 
galaxy (see Caranicolas & Zotos 2009b, 201 la, 201 lb). The 
first real surprise in the study of dark matter lay in the out- 
ermost parts of galaxies known as galaxy halos. Here, there 



is negligible luminosity, yet there are occasional orbiting gas 
clouds, which allow us to measure rotation velocities and dis- 
tances. The rotation velocity was found not to decrease with 
increasing distance from the galactic center, implying that the 
mass distribution of the galaxy cannot be considered like the 
light distribution. The mass must continue to increase since 
the rotation velocity satisfies 



where M is the mass within radius r\ we infer that M in- 
creases proportionally to r. This rise appears to stop at about 
50 kpc, where halos appear to be truncated. We infer that 
the mass-luminosity ratio of the galaxy (M/L), including its 
disk halo, is about 5 times larger than estimated for the lu- 
minous inner region or equal to about 50. The rotation curve 
shown in Figure 3, corresponding to our disk galaxy, seems 
to maintain high values of the angular velocity for large 
distances from the galactic center. Thus, we may say that our 
dynamical model can interpret the presence of dark matter. 

Forty years ago, galactic activity and interactions be- 
tween galaxies were viewed as unusual and rare. Nowa- 
days, they seem to be segments in the life of many galaxies. 
From the astrophysical point of view, in the present work, 
we have tried to connect galactic activity and galactic inter- 
actions with the nature of orbits (regular or chaotic) and also 
with the behavior of the velocities of stars in the primary 
galaxy. We consider the outcomes of the present research 
to be an initial effort, in order to explore the dynamical struc- 
ture of the 3D binary stellar system in more detail. As results 
are positive, further investigation will be initiated to study 
all the available phase spaces, including orbital eccentricity 
(elliptic orbits) of the small companion galaxy and its incli- 
nations to the primary galaxy. Moreover, we shall try to use 
the outcomes obtained from this initial and simple dynamical 
model to conduct computer yV-body simulations in a binary 
system of interacting galaxies (a primary galaxy with a satel- 
lite companion), in order to reveal changes in their orbital 
properties through merger processes and tidal effects, which 
are obviously out of the scope of this study. 
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